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Abstract 

In recent years, cell population models have become increasingly common. In contrast to classic single cell models, 
population models allow for the study of cell-to-cell variability, a crucial phenomenon in most populations of 
primary cells, cancer cells, and stem cells. Unfortunately, tools for in-depth analysis of population models are still 
missing. This problem originates from the complexity of population models. Particularly important are methods to 
determine the source of heterogeneity (e.g., genetics or epigenetic differences) and to select potential (bio-) 
markers. We propose an analysis based on visual analytics to tackle this problem. Our approach combines parallel- 
coordinates plots, used for a visual assessment of the high-dimensional dependencies, and nonlinear support 
vector machines, for the quantification of effects. The method can be employed to study qualitative and 
quantitative differences among cells. To illustrate the different components, we perform a case study using the 
proapoptotic signal transduction pathway involved in cellular apoptosis. 



1 Introduction 

Cell populations are heterogeneous in terms of, e.g, cell 
age, cell cycle state, and protein abundance [1,2]. This 
heterogeneity is ubiquitous, even in clonal population, 
and influences cell fate decisions [2,3], such as cell 
death/proliferation [4-7]. Thus, to ultimately understand 
and control the behavior of populations, the key sources 
of cell-to-cell variability have to be unraveled. Unfortu- 
nately, this is challenging due to experimental con- 
straints. Most experimental systems and measurement 
devices only allow for the simultaneous assessment of a 
few cellular properties on a single cell basis. This prohi- 
bits the purely experimental analysis of processes which 
depend on many different cellular properties. Spencer et 
al. [5] have shown that the experimental limitations can 
be overcome partially using mathematical models. 

To mathematically describe heterogeneous popula- 
tions, agent-based models are used most frequently. 
Each agent provides a mechanistic description of the 
signal transduction within individual cells and thus of 
its behavior. In such a framework, variability can be 
modeled by either stochastic [8-10] or deterministic 
[4,5,11] differences among individual cells. The source 
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of the former is the stochasticity of biochemical reac- 
tions, while the latter may arise from genetic and epige- 
netic differences, environmental heterogeneity, or slow 
dynamic processes (such as the cell cycle). 

We focus on the deterministic differences among cells 
— also called extrinsic factors [12] — in populations of 
non-interacting cells. Those differences are commonly 
modeled by differential parameter values and initial con- 
ditions [5,13]. Several methods exist to infer the distri- 
bution of parameters and initial conditions from 
experimental data [13-15] and to obtain quantitative, 
mechanistic models for cell populations. Unfortunately, 
the resulting agent-based models are in general highly 
complex. This complexity prevents the analysis of these 
models using common tools for dynamical systems [16], 
such as sensitivity and bifurcation analysis. To the best 
of our knowledge, for models of heterogeneous cell 
populations, no structured analysis approach is available. 
To study population models and to facilitate a model- 
driven analysis of the heterogeneity, highly flexible 
methods are required which do not rely on an analytical 
analysis. 

In this work, we propose two methods to fill this gap 
and to facilitate the analysis of population models. 
These methods — parallel-coordinates plots [17] and 
support vector (SV) machines [18-20] — are tools widely 
used for the analysis of high-dimensional datasets. We 
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outline how these tools can also be used to analyze 
complex models of heterogeneous cell populations, par- 
ticularly addressing the question: "Which parameters 
cause the heterogeneity of the population's response?". 
Thereby, we extend our previous work [21] and consider 
qualitative heterogeneity among cells, in the context of 
cell fate decisions, as well as quantitative heterogeneity, 
such as the delay of a decision process. 

We show that parallel-coordinates plots provide an 
easy tool to obtain a qualitative understanding of the 
system, whereas SV machines allow for assessing the 
performance of marker combinations quantitatively. 
Good markers are thereby defined as single cell para- 
meters that facilitate a good prediction of the cell fate 
decision or the quantitative property under considera- 
tion of the individual cell. Furthermore, we show how 
the combination of parallel-coordinates plots and SV 
machines enables an in-depth analysis of complex mod- 
els using exploration techniques. 

The article is structured as follows: In the section 
"Methods", the considered system class and problem are 
described in mathematical terms, the general idea is dis- 
cussed, and the application of parallel-coordinates plots 
and SV machines is outlined. In the section "Results", 
we provide an exemplary application of our method to a 
model of the caspase cascade. The article is summarized 
in the section "Discussion". 

2 Methodology 

2.1 Models for heterogeneous cell populations and 

decision processes 

2.1.1 Mechanistic population model 

In this article, population dynamics are described using 
an ensemble [5,13] of cells (agents). This yields the 
agent-based population model: 

£ pO p = {£(0»)|i = {l N], 0«~©(0)}, 

in which the superscript (i) specifies individual cells 
within the population, N e N denotes the size of the 
cell ensemble and £(t?' !) ) is the model of the 2-th cell. 
The sine le cell model T.(9 {i) ) may belong to the class of 
Markov jump processes [15], stochastic differential 
equations [14], or ordinary differential equations [13]. 
Since in this study we are mainly interested in signal 
transduction and decision making, we consider ordinary 
differential equation models. Each individual cell of E p0 p 
is described by 

E(0«] :x® = f{x {l \ 9% x {i] {0)=x Q {9^), 

with state vector x^(t) e R" and parameter vector 
0(0 e Rf . The vector field / : R» x R? -> R" describing 



the cell dynamics is locally Lipschitz and the mapping 
Xo : R? —> R" is continuously differentiable. The para- 
meters 6^ may be kinetic constants, such as synthesis, 
degradation, or reaction rates. 

Heterogeneity among cells of the ensemble is modeled 
by differential parameter values (r and initial conditions 
x 0 {8^) among individual cells. The density of parameters 
9 is given by a probability density function 
© : -> R + . Thus, the probability of observing 9 (i) e D. 
is 

Prob(6»W g Q) = j &(9)d9. 

Q 

This modeling framework is highly flexible and has 
been proven to be very useful, especially if fast signal 
transduction processes, such as cellular apoptosis, are 
investigated. For a more detailed introduction, we refer 
to the work of Spencer et al. [5] and Hasenauer et al. 
[14]. The properties of such populations of single cells 
have been studied by Spencer et al. [5], while Hasenauer 
et al. [14] have derived a partial differential equation 
model for the resulting population dynamics. 
2.1.2 Qualitative and quantitative properties of the single 
cell response 

Given the mathematical models introduced above, we 
study qualitative and quantitative properties of the single 
cell responses. Qualitative properties are defined as the 
outcome of a discrete decision processes, e.g., whether 
the state of a bistable system converges to one or the 
other stable steady state, or whether a certain concen- 
tration threshold is reached. In contrast, quantitative 
properties allow the assessment of small differences 
among cells, such as the time point when a particular 
threshold is exceeded. 

To define single cell properties given the single cell 
trajectory x {l \-), the functionals F^: i 1 -> R and F$: I 1 — > 
{-1, +1} are introduced. The functional is used to 
evaluate the quantitative property cfr = F^x (•)) £ R, 
while Fs determines the qualitative property <r'' = Fg{x^ 

(•)) e (-1. +!}• 

To exemplify the functionals, we consider a process in 
which threshold exceeding and its timing are of interest. 
Such processes are important, for example, in apoptotic 
signaling [5] and cell cycle progression [22,23], and 
allow for two outcomes. Either the concentration of a 

molecule Xj within the 2-th cell exceeds the threshold 

Xi t h, <5 = +1, or it does not, 8^' = -1. This yields the 
decision functional 

F S {^{-)):=\ +1 ^H 0 ^** (1) 
— 1 otherwise. 
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For the subgroup of cells exceeding the threshold, the 
time of threshold exceeding is defined by the second 
functional 

i>(x w (-)) := argmin jxfft) > x jlth \ , (2) 

and may be employed to achieve a quantitative 
understanding. 

Note that the response # (•) of a cell merely depends 
on the cell's parameters (9 W , as the single cell model is 
deterministic. Therefore, the quantitative and qualitative 
properties of a single cell can be viewed as a function of 
the parameters, 0 (o = <j>{6 (i) ) and 8 (i) = 8{& i] ). Differences 
in the parameters — as they arise between different cells 
— may hence influence 8 and (f> , which determine 
cell fate decision and qualitative properties of the cells. 
2.1.3 Response markers 

To understand the heterogeneity within the population 
response S p0 p, it is necessary to assess the dependency 
of d and (j> on the individual parameters 0j. In parti- 
cular, the question arises which subset 6 m of para- 
meters, 

&m ■= [Bm v dmf, with m C {1 {ft, 

is responsible for which aspect of the population het- 
erogeneity. Mathematically, m is an index set and, e.g., 
for m = [2, 4] T only 0 m = [6 2 , 6> 4 ] T is considered. The 
question of the relative importance of different para- 
meters directly relates to the common problem of bio- 
marker selection for stem cells and tumor cells, which is 
experimentally challenging. 

If there exists a subset 6> m of the parameters 8 which 
allows for the reliable prediction of the response, not all 
sources for heterogeneity have to be assessed but only 
those associated to f? m . This enables a focusing of the 
model development, as well as the reduction of the 
experimental effort. 

2.2 Analysis of population models using data analysis 
tools 

In this contribution, we illustrate the application of par- 
allel-coordinates plots and support vector machines for 
the study of parameter dependencies and the selection 
of markers m. Parallel-coordinates plots and SV 
machines are well-known, but almost exclusively applied 
to study high-dimensional sets of measurement data. To 
exploit the methods for the analysis of simulation mod- 
els, at first the cell ensemble is simulated for N » 1. 
This yields many pairs of parameters and trajectories, 

(>>, i=l, N, 



which are then used to obtain samples of quantitative, 

S 9 = ((V'V) (o™, <pW)}, withpW = F p (rM(.)), 

and qualitative 

S s = {(0M,« (1) ),."'( flC "°' a( ' 0 )}' withSffl =Pj(x«(0), 

cell properties of interest. These samples contain 
information about the dependency of (f> and 8 on the 
parameters 6, being analyzed in the following. To study 
the high-dimensional mappings 8 = 8(9) and 4> = 4>(9), 
parallel-coordinates plots will be employed. For the 
quantitative assessment of particular marker combina- 
tions SV machines will be applied. By combining both 
approaches it is possible to quickly gain an overview of 
important interrelations and quantify those. 
2.2. 1 Combining parallel-coordinates plots and SV machines 
to a visual analytics system 

The proposed simulation data-based analysis approach 
circumvents an analytical analysis of the system equa- 
tions, which would be time consuming and could only 
be carried out by experts. However, the simulation data- 
based approach creates the need for analyzing the large, 
high-dimensional datasets, Ss and S v . 

The analysis of such datasets often relies on a reduc- 
tion of complexity while preserving the important infor- 
mation. Visualization can help in such a situation to 
determine the important parameters and to avoid infor- 
mation loss. In this work, parallel-coordinates plots are 
used to gain insight into the high-dimensional depen- 
dencies and to find interesting dimensions. In this parti- 
cular setting, interesting dimensions are those that 
clearly separate a given set of classes and thus are good 
candidates for the selection of potential markers m. In a 
second step, the potential markers m are used to train a 
SV machine. These SV machines allow for a quantitative 
evaluation of the marker quality. While SV machines are 
also helpful on their own, checking all possible combi- 
nations of markers would result in a combinatorial 
explosion. By combining SV machines and parallel-coor- 
dinates plots, the number of necessary SV machine eva- 
luations can be decreased substantially, resulting in a 
tremendously reduced computational complexity. The 
overall workflow of the analysis illustrated in Figure 1. 

Besides an improved understanding of the model, 
results obtained during the analysis can be used to 
adapt the population model or to select additional 
experiments. This proposed framework, integrating 
interactive visualization with automated methods while 
allowing for a feedback to the actual system/model, thus 
incorporates important aspects of visual analytics [24]. 
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Experiments / experimental data 
(e.g., flow cytometry, microscopy, . . .) 



Modeling and parameter estimation 

(see, e.g., Koeppl et al. (2011) and Hasenauer et al. (2011a, b)) 
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Figure 1 Visual analytics for marker selection. Illustration of the overall workflow of (1) experiments and the collection of measurement data, 
(2) modeling and parameter estimation, and (3) model analysis. Our visual analytics framework for marker selection in models of heterogeneous 
cell populations is shown in detail. Based on simulation data Ss and S v generated from the model, a visual analysis is performed using 
parallel-coordinates plots. Employing this visualization, insight can be gained into the dependencies of the considered properties on the 
parameters. Additionally, a potential marker combination m can be selected, here, for instance, m = 1 or m = 4, which allow for a separation of 
the classes. For the quantitative assessment of the informativeness of 9 m , e.g., SV classification or SV regression, may be employed. 



2.3 Parallel-coordinates for the analysis of high- 
dimensional data 

Parallel-coordinates [17] are a popular visualization 
technique for high-dimensional data. A parallel-coordi- 
nates plot is constructed by placing axes in parallel, as 
illustrated in Figure 2. A single pair of adjacent axes 



represents a 2-D projection of the data, where a point of 
the corresponding Cartesian coordinates is mapped to a 
line in parallel-coordinates, and vice versa. Due to this 
point-line duality, the same patterns emerge in a paral- 
lel-coordinates plot as in the dual Cartesian coordinates. 
However, adding more axes not only allows to visualize 
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xi x 2 x 3 



Figure 2 Parallel-coordinates plots represent multi-dimensional data as polylines crossing parallel axes A point in Cartesian coordinates 
is mapped to a line in parallel-coordinates. As more axes are added, a line can be traced over all dimensions, which greatly facilitates the 
perception of multi-dimensional data. 



a set of pairwise relations, but also supports the viewer 
in tracing lines over all dimensions. As a result, multi- 
dimensional outliers and clusters can be visualized 
together with 2-D relations and the distribution of 
values for single dimensions. 

As an iV-dimensional data point is represented by a 
polyline intersecting axes at the respective values, paral- 
lel-coordinates greatly suffer from overplotting if many 
lines have to be drawn. In the resulting clutter of lines, 
interesting patterns might be hidden from the user. 
Exploiting the point-line duality, similar clutter-reducing 
approaches as for Cartesian coordinates can be used, 
where a popular technique is to estimate the density of 
points (lines) and to render points (lines) transparently 
with blending enabled. Other approaches compute a 
continuous density [25] or estimate the overall density 
using density estimation techniques [26,27]. In this 
work, both alpha and additive blending is used to visua- 
lize the parameter distribution in the different classes 
((jy® = 1 and </> w = -1), enabling a qualitative analysis of 
their multi-dimensional shape. An example of this alpha 
blending is shown in the section "Results". 

For the analysis of a continuous variable, colormaps 
can be applied to the axis representing the dependent 
variable (fr'\ Then, every polyline is rendered using a 
color according to <f> , such that its value can be visually 
determined over the whole plot. The overall distribution 
of colors can then be used in conjunction with the 
shape of lines to analyze the dependency of independent 
variables from the dependent. Again, overplotting can 
become an issue for large datasets, such that a separa- 
tion in few classes and a separate visualization of those 



might be more informative (see example in section 
"Results"). 

2.4 SV machines for the quantification of marker 
performance 

Given a basic understanding of the importance of the 
parameters and a potential marker combination 6 m , a 
quantitative assessment of the predictive power of £> m is 
desirable. To achieve this, the samples S$ and S v are 
analyzed employing nonlinear SV classification and non- 
linear SV regression, respectively. SV classification allows 
for the study of decision processes, while SV regression 
enables the analysis of quantitative system properties. 

The performance of SV machines — which might be 
interpreted as data-based predictors — provides a mea- 
sure for the quality of the marker combination # m . If a 
SV machine using only (9 m provides good predictions for 
a decision process which depends on 9, then this means 
that 9 m carries the most important information. This 
will be discussed in more detail in the following. 
2.4. 1 SV classification 

The goal of the SV classification is to predict the dis- 
crete property (5 W given Q W . Thus, the nonlinear map- 
ping 8 = 8{ff) is approximated by the lower-dimensional 
nonlinear mapping 8 = S(0 m ) ■ To calculate the SV clas- 
sifier, a two step procedure is applied, as illustrated in 
Figure 3. First, a mapping <j> : R r — >■ R r *— also called 
kernel — is constructed that transforms the input space 
into a feature space of higher dimension (r* >r). Second, 
a linear separation of the data is performed in feature 
space [20]. Therefore, the optimization problem 
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Figure 3 Support vector machine for subgroup classification 

Visualization of the SV machine approach for separating cells with 
<5 W = +1 (+) and <S (,) = -1 (°). Left: distributed data in the input space. 
Middle: sample transformed in the feature space which allows for 
better separation. Right: separation result for separating hyperplane 
with normal vector w (— >). As a perfect separation is in general not 
possible, misclassifications (o) exist. 



minimize — w w + 

w,b,( 2 



subject to <5<'> (w T 4>(8£) + b\ > 1 - fj, i = 1, . . . ,N, 

?j > 0, 



(3) 



i = 1, 



is solved, in which w and b denote the normal vector 
of the separating hyperplane and its offset, respectively. 
The objective function combines a misclassification pen- 
alty, Yli=i ti ' anc ^ a mar g m maximization, jW T w ■ The 
weighting of the different terms can be influenced via C. 
The constraints are that all data points (S>(9$,8®) are 
correctly classified within a certain error margin iff,. 

Given the solution of (3), a predictor (SV classifier) for 
the decision process S = 3(9) is 



5(0 m ) = sign {w T <t>(9 m ) + b). 



(4) 



Assuming that the training set Sg is large, the predic- 
tive power of this predictor will be high — meaning that 

8(6$) = 8(0®) for most 6 {i) ~ 0(0) — if and only if the 
selected markers 6 m are informative. This allows the 
quantitative assessment of the informativeness of the 
markers 6 m using the SV classifier. 

Therefore, a second sample S' s is computed which 
was not used to train the SV classifier, avoiding overfit- 
ting. For this sample, the predictor 5(0 = 8(9$) is eval- 
uated. These results are used to calculate the percentage 
of true positive classifications TP (5(0 = i /\ j(0 = 1) 
and false positive classifications FP (5(0 = o A 8® = 1) 
achieved by the SV classifier. TP and FP provide infor- 
mation about the predictability of the outcome for 0 
using solely q$ . Thus, the marker quality can be 
assessed via TP and FP. If a low-dimensional m exists 
that provides TP ~ 1 and FP » 0, the parameters 8 m 
dominate the decision process and are good markers. 



For a quantification of this effect, the classification per- 
formance can be analyzed in receiver-operating charac- 
teristic (ROC) space [28]. 
2.4.2 SV regression 

Similar to the assessment of the predictive power of 
marker combinations for qualitative decisions, also 
quantitative properties may be analyzed. Therefore, we 
employ SV regression which allows us to compute a 
data-based predictor 



<p(9 m ) = w T <t>(9 m ) + b, 



(5) 



for the quantitative property <f> = (f) (6). To compute 
the nonlinear predictor, a kernel <J> ; [R r — > [R r * [29] is 
chosen and an optimization criterion selected. In this 
work, we use an e-insensitive loss function [30], mean- 
ing that residuals ^(0 _ <p(9$) with an absolute value 
below £ are not penalized while larger residuals are 
penalized linearly. This loss function is frequently used 
in the literature (see, e.g., [20,30]) and results for the 
sample 5^ in the optimization problem: 



N 



minimize —ww 

w,b££* 2 



subject to <p 



(0 



+ C + 

i=i 

w T t>(0$) 



-b> e + | i; 
-<p® +w T <t>(9$) 
+b>s + £*, 



1,...,N, 



(6) 



i = 1, 
i = 1, 



,N, 
,N. 



Aside from the penalization of prediction error, 
zL^ife + £*), flatness and a unique solution is ensured 

using jW J w. The trade-off between those two is deter- 
mined by the constant C > 0. 

The optimal solution of (6) for w and b provides the 
optimal predictor (5) with respect to the loss function 

and kernel. This predictor ^(0 = (p(9$) is applied to a 
second sample S' v to compute r/>(0, a prediction for <jr l \ 
Employing <fr 1 ' and r/>(0 the marker combination m 
might be evaluated based on the relative prediction 



errors, e, 



(0 



Using g(0, the prediction 



powers of different marker combinations can be 
assessed and compared using, e.g., the mean error 

i 6m • ^ tne mean prediction error achieved by a 
marker combination is small, the parameters ff$ carry 
most of the information about (fr*\ and hence are suita- 
ble markers. In some situations, the information about 
the mean prediction error may be complemented by 
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detailed information about the error statistics, {eS}^- 
These statistics may be visualized using, for instance, 
box plots or histograms, and provide additional insight, 
e.g., in the structure of the error (short- vs. long-tailed 
distributions) and the potential causes. 

Note that the performance and predictive power of SV 
machines strongly depend on the available training set. 
For the analysis performed, we ensured that the training 
sets are large enough and that a further increase in its 
size does not result in a significant improvement of the 
predictors. This is, in most situations where SV 
machines and SV regressions are used, impossible for 
data analysis, as the measurement devices are limited. 
However, in this work we study the problem of model 
analysis. The size of the dataset can be increased arbi- 
trarily by repeated simulation of the model. Besides the 
size of the dataset, the parameters of the SV classifica- 
tion and SV regression are tuned to allow for a fair 
comparison between the marker combinations. With 
this and the existence of sophisticated SV machine tool- 
boxes (e.g., LIBSVM [31]), the observed difference 
between marker combinations can be assumed to be 
due to the predictive power of the markers. 

Summing up, SV machines allow for the derivation of 
predictors for qualitative and quantitative properties. 
These predictors can be used to assess the information 
content of a subset m of the parameters about the 
respective properties, thereby facilitating the assessment 
of a quantitative evaluation of the predictive power of 
0 m . For further details about SV machines we refer to 
[18-20,30,31] and references therein. 

3 Results 

3.1 Model for heterogeneous cancer cell population 

To illustrate the proposed visual analytics framework, a 
model of the proapoptotic signaling is analyzed. Proa- 
poptotic signaling is involved in the process of apoptosis 
[32-34], also called programmed cell death. Apoptosis is 
an important physiological process to remove infected, 
malfunctioning, or no longer needed cells from a multi- 
cellular organism. The apoptotic signaling pathways 
converge at the caspase cascade [32], where initiator 
caspases (e.g., caspase 8) and effector caspases (e.g., cas- 
pase 3) are activated. If the activity of effector caspases 
exceeds a certain threshold, apoptosis is induced. 

A variety of single cell and cell population models 
have been proposed to describe cellular apoptosis (see, 
e.g., [4-6,34-40] and references therein). In this study, 
we consider the mathematical model of the signal trans- 
duction which is depicted in Figure 4. This single cell 
model [35] is among the most cited ones. For details 
about the model, we refer to the original publication 
[35]. As the process of apoptosis induction is known to 



stimulus 

(e.g. TNF or TRAIL) 



C8 



" 0 
IAP 





C3 



cell death 

Figure 4 Caspase cascade. Illustration of proapoptotic signaling 
pathway [35]. Normal arrows refer to conversion reactions, dashed 
arrows indicate enzymatic activity, and thick arrows illustrate inputs 
and outputs of the system. 



be heterogeneous, we extend the single cell model [35] 
by accounting for cell-to-cell variability. This is achieved 
by introducing differences in parameter values and 
initial conditions: 

♦ From flow cytometric experiments, it is known 
that the amount of caspase 8 (C8), caspase 3 (C3), 
caspases 8- and 10-associated RING protein (CARP), 
and inhibitor of apoptosis protein (IAP) is different 
among individual cells. The differences are modeled 
by differences in synthesis rates (k_ & , /c 9 , k. w , and k. 
12 ) among individual cells. The distribution of /e 8 , k. 
9, /cio, and k_i 2 within the population is modeled as 
log-normal distribution, with mean as published by 
Eissing et al. [35] and a coefficient of variation of 0.4 
(own unpublished data). The initial conditions of C8, 
C3, CARP, and IAP are set to their steady state 
values. 

• Similar to the original publication [35], the activa- 
tion of the caspase cascade is modeled by a non- 
zero initial condition of active caspase 8, C8a(0). In 
the population, C8a(0) is log-normally distributed 
with a median of 4,000 molecules per cells and a 
coefficient of variation of 0.4. The variation of C8a 
(0) accounts for variability up-stream of the caspase 
cascade. 

The binding affinities and kinetic rates are the same 
for all cells. For the numerical values, we refer to the 
article of Eissing et al. [35]. 

Given this model of the heterogeneous cell population, 
we analyzed (i) how the decision whether or not a cell 
undergoes apoptosis during the first 12 hours and (ii) 
how the time of cell death T d is influenced by the cell's 
parameters 9 = [C8a(0), k.g, ks, k-w, £-12] 1 • This yields 
two variables of interest: S (= +1 => cell survived; = -1 => 
cell died) providing the outcome of the decision process; 
and 4> (= T d ) providing the time of apoptosis commit- 
ment. As indicator for apoptosis, the amount of active 
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caspase 3 (C3a) is used. If more than 5,000 copies of 
C3a are present in a cell, this cell is assumed to undergo 
apoptosis within 10 minutes, defining the time of cell 
death T d . The functionals associated to the considered 8 
and 4> are similar to (1) and (2), respectively. In the 
remainder, we search for a lower-dimensional subset of 
the parameters 9 which provide good markers for cell 
death and survival as well as the time of cell death. 

3.2 Parallel-coordinates plot establishes importance of C3 
and IAP concentration for cell fate decision 

To study the life-death-decision, a sample S$ with 
100,000 members is visualized in parallel-coordinates 
(Figure 5). As only two classes (dead and alive) are con- 
sidered, alpha blending can be used to visualize the den- 
sity of each class as well as the density at the 
overlapping regions, where the transparent red color, 
representing dead cells, and the transparent blue color, 
representing living cells, are blended wit a = 0.03. Using 
this coloring, high-density regions appear more satu- 
rated for the individual classes and darker at their 
overlap. 

From Figure 5, it is apparent that the second and 
fourth parameters (f? m = [k. s , /c_io] T ) provide a reason- 
able separation between the classes (red = dead, blue = 
alive). Most of the surviving cells have high values of k.g 
and low values of /ci 0 , which corresponds to a high IAP 
expression and a low C3 expression, respectively. 
Although the other parameters also influence the pro- 
cess, their influence seems to be minor. 

3.3 SV classification proves that C3 and IAP expression 
are the best markers for the cell fate decision 

Given the results of the visual analysis, we consider d m 
= d m = /c_io, as well as 8 m = [k_ s , k. w ] T and compute 
the classification quality using SV machines (for details 
see "Methods"). As can be seen in Figure 6A, the 



predictive power of the individual parameters is limited 
(0 m = L s : TP = 0.73, FP = 0.38; 9 m = L 10 : TP = 0.74, 
FP = 0.29), while both markers together yield a reason- 
able classification performance (TP = 0.77, FP = 0.13). 
The corresponding ROC curve is depicted in Figure 6C 
and the visualization of TP and FP is provided in Figure 
6D. For comparison, the alternative combinations of two 
markers are evaluated in terms of the area under the 
ROC curve (Table 1) and the TP/FP (Figure 6C). 

The markers 6 m = k. s and (9 m = /c 10 outperform all 
other single markers and marker pairs. In addition, the 
marker vector 6 m = [k. s , k. l0 ] T outperforms all other 
combinations in terms of the area under the ROC 
curve. Some other combinations result in more than 
50% false positive classifications (see Figure 6B). Of 
course, extending the marker vector, e.g., by adding k^ 2 > 
results in further improvement. 

3.4 Parallel-coordinates plots show a complex 
dependency of the time of death on the parameters 

After the analysis of the decision process, we study the 
dependency of time of cell death T d on the parameters. 
The time of cell death T d is a quantitative property and 
can take any positive value, therefore an alternative 
visualization has to be used. One approach would be to 
use a different color for each line in parallel-coordinates, 

depending on jjp. Unfortunately, this approach suffers 

from heavy overplotting, which is why the data was split 
into three classes and separate plots were created for 
each class. 

Figure 7A-C visualize the parameter distribution in 
different percentile intervals for T d . A comparison of 
Figure 6 A, visualizing the cells that die early (0 to 10th 
percentile), and Figure 7C, depicting the cells that die 
late (90 to 100th percentile), unravels o sets in all para- 
meter dimensions. The differences are particularly pro- 
minent for C8a(0), k_io, and k_ 12 , showing that the 




C8a(0) ks fc-9 fe-io fc-12 

Figure 5 Separation of subgroups in parallel-coordinates. Parallel-coordinates density plot in which each polyline represents the parameter 
of a single cell, 6> ( ". The color of a polyline encodes whether the cell survived (blue line) or died (red line). In order to emphasize dense regions, 
alpha blending with a = 0.03 was used for all lines. The parameters k. a and k. ]Q show the best separation of colors and hence correspond to 
potential markers. 
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IAP synthesis, fc_8 



initial active C8, C8a(0) 





C8a(0) and fc_( 
C8a(D) and fc_( 
C8a(0) and k _, 
C8a(0| and fc_, 



False positive False positive 

Figure 6 Performance of different marker combinations. Evaluation of classification performance using different marker combinations. (A) 
and (B) illustrate the classification obtained using to different two marker combinations. The prediction of the classifier (trained using 10.000 
simulations) is provided as background color (blue square = alive; red square = dead) together with a test sample ( x = alive; ° = dead), which 
has not been used for training. (A) Classification employing C3 synthesis, /C| 0 , and IAP synthesis, k. 8l as markers. For the classification of cell 
survival: TP = 0.77, FP = 0.13. (B) Classification employing initial amount of C8a, C8a(0), and C8 synthesis, /;_ 8 , as markers. For the classification of 
cell survival: TP = 0.68, FP = 0.53. (C) ROC curve obtained when using C3 synthesis, /c 10 , and IAP synthesis, fc 8 as marker. (D) The performance of 
all individual markers, all marker pairs, and the best marker triplet in ROC space. Note that an optimal classifier would be in the upper left corner. 



abundance of C3 also plays an important role in deter- 
mining whether cells die early or late. Unfortunately, a 
closer look at Figure 7 also reveals that the parameter 
distributions associated to cells that undergo apoptosis 
at early, intermediate, and late time points strongly 
overlap in parallel-coordinates. This indicates that T d 
may depends on all parameters. Therefore, a reliable 
prediction of Td using only a few parameters might be 
infeasible. 



Table 1 Area under the ROC curve for different marker 
combinations 





C8a(0) 




k-9 


'f-io 




C8a(0) 


0.569 


0.747 


0.626 


0.808 


0.690 


k. a 


0.747 


0.736 


0.760 


0.898 


0.800 


k-9 


0.626 


0.760 


0.603 


0.822 


0.709 




0.808 


0.898 


0.822 


0.795 


0.858 




0.690 


0.800 


0.709 


0.858 


0.676 



A common performance measure for classifiers is the area under the ROC 
curve. The worst and best achievable performance is 0.5 and 1, respectively. 
The evaluation of the area under the ROC curve verifies that if only one 
marker can be used, /V_ 10 (area = 0.795) and k. 8 (area = 0.736) are the best 
choice (on-diagonal bold numbers). In case of two markers, the combination 
of /c_ 8 - /c„i 0 allows for the best classification (area = 0.898) (off-diagonal bold 
numbers). For the three markers /c_ 8 , /f_i 0 , and /c_ 12 , the area under the ROC 
curve is 0.966. 



3.5 SV regression reveals ubiquitous importance of IAP 
an C3 expression levels 

To quantify the predictive power of different marker com- 
binations with respect to T d , we employ the SV regression 
based approach introduced in "Methods". As a perfor- 



4 ] - 



mance measure, the relative prediction error 



their fy is the prediction of the SV machine. Details on 

the implementation may be found in "Methods". 

At first, we study the potential combinations of two 
markers proposed by the parallel-coordinates plots: k.io 
and /c_i 2 ; C8a(0) and £_ 10 ; and C8a(0) and k. u - Out of 
those, the best performance with a median prediction 
error of 40% is achieved by C8a(0) and /c 12 , which also 
outperforms all other combinations of two markers. Inter- 
estingly, all marker combinations achieve a median predic- 
tion error between 40 and 50%, as shown in Figure 8. This 
illustrates two things: On the one hand, markers allowing 
for a distinction between early and late dying cells do not 
necessarily enable a good prediction of the death time T d , 
as here also the cells dying in an intermediate interval 
dominate the statistic. On the other hand, this quantifica- 
tion proves that even the best combination of two markers 
provides only very limited predictive power. Thus, unlike 
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c 




C8a(0) fc_ 8 fc_ 9 fc_io fc_i 2 

Figure 7 Dependency of T<j on the single cell parameters. Visualization of three subsets of the sample iS p using parallel-coordinates. The 
subplots depict the parameter vectors of the cells having value: (A) below the 10th percentile; (B) between 45th percentile and the 55th 

percentile; and (C) above the 90th percentile of the T d values. The thin black lines are elements of the sample and the thicker line is the mean. 
A comparison of the subplots (A) and (C) shows that cells dying early and cells dying late mainly differ in the parameters C8a(0), /c 10 , and /c 12 
but also the other parameters show small offsets in one or the other way. The joint consideration of all subplots reveals that in parallel- 
coordinates the parameter regions, associate to different times of cell death, overlap. This indicates that several parameters will determine TV. 



the decision which predominantly depends an C3 and IAP 
expression, the time of cell death is highly sensitive to 
changes in all parameters. 

4 Conclusion 

4.1 Visual analytics enable an in-depth analysis of 
complex population models 

In this article, a novel explorative approach has been 
presented to determine markers for decision processes 



in heterogeneous populations. It has been shown that 
methods used for data analysis can also be employed to 
gain insight into complex models, where common analy- 
tical methods seem to reach their limits. Especially, the 
potential of parallel-coordinates plots and support vec- 
tor machines has been illustrated. While the first allows 
for the study of large, high-dimensional datasets and the 
selection of potential markers, the latter can provide a 
quantitative assessment of their predictive power. Using 
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Figure 8 SV regression analysis of Tj. Statistical analysis of the predictive power of SV machines using different markers m. Shown are box 
plots for the relative prediction error: median (blue line), interquartile range (light blue area), and 10th to 90th percentile range (green line). The 
combination of two markers which results in the smallest median prediction error is (C8a(0) and /c 12 ), but all over combinations achieve 
prediction errors which are alike. The last column shows that information about all parameter enables extremely precise predictions using SV 
machines. 



both methods, the source of qualitative and quantitative 
cell-to-cell variability may be unraveled. 

This article provides a case study evaluating the 
potential of combining visualization and automated 
methods for the assessment of complex system models. 
The considered system class is only one example and 
the proposed framework can be generalized easily to 
other systems and questions. 

4.2 Analysis of heterogeneous cell population allows for 
novel insight 

We have illustrated the proposed visual analytics approach 
by analyzing a cell population model for proapoptotic sig- 
naling, which plays an essential role in programmed cell 
death. We have studied the cell fate decision as well as the 
time of cell death. These properties were analyzed before 
(see, e.g., [5]) in a purely qualitative way and without the 
tools proposed in this work. 

Our study shows that parallel-coordinates plots are a 
proper tools to determine potential markers. The 



predictive power of these markers can then be quantified 
using SV machines. In this study, the markers we found 
agree well with those found in the literature. In particu- 
lar, the important role of IAP — also called XIAP — for cell 
death commitment is outlined in several publications 
[39,41]. While C3 abundance is known to be important 
[39], our analysis suggests that the amount of available 
C3 could be even more important than expected. 

In addition, our analysis indicates that, under normal 
conditions, the time of cell death strongly depends on 
all parameters, which has been hypothesized earlier [5]. 
Only under altered conditions, e.g., a strongly increased 
initial amount of C8a(0), some parameters become more 
important than others (results not shown). This is again 
in agreement with the results of Spencer et al. [5]. 
Furthermore, this finding of a varying importance of 
parameters depending on the experimental setup, pro- 
vides hints for possible future experiments. Thus, our 
visual analytics approach we propose also provides help- 
ful feedback for model validation and development. 
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4.3 Outlook 

In this work, we have proposed a method to determine 
decision markers for given models. However, all model 
possess uncertainties, rendering an uncertainty-aware 
analysis crucial. Therefore, a workflow including model 
development, parameter estimation, uncertainty analysis, 
and marker prediction has to be established. This 
requires improved modeling and parameter estimation 
tools, as well as methods to evaluate the uncertainty of 
the marker prediction, arising from model uncertainties. 

Given such a workflow, beyond the analysis of models, 
our analysis tools might also be used to guide the search 
for biomarkers. This is possible as our methods allows 
for the assessment of the importance of any parameters 
which are different among cells of the population. 
Among others, the importance of common biomarkes, 
e.g., expression levels and transcription factor/protein 
abundance, may be determined based on a model of the 
population. This is much in the same way as the target 
selection using sensitivity analysis of single cell models 
based on ordinary differential equations (see, e.g., [42]). 
However, the marker selection requires population mod- 
els, as differences between cells have to be considered, 
and is therefore more challenging. 

Methods 

Software 

The model of the heterogeneous cell population was 
implemented in MAT LAB using the SBtoolbox2 [43]. 
For the SV classification and the SV regression, the 
LIBSVM toolbox for MATLAB is employed [31]. The 
visualization software for parallel-coordinates was imple- 
mented in C++ using the Qt library Version 4.8.0 and 
OpenGL. 

Numerics 

For the SV classification and SV regression, we 
employed as kernels radial basis function with y = 0.25. 
The SV regression parameter which defines the interval 
of insensitivity was set to 0.01. All remaining parameters 
are set to the default constants, see LIBSVM manual. To 
improve the performance of the SV machines, we 
applied a log-transformation to the parameters 0. 
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